Data analysis
Vocal interactions in Thyroptera tricolor
Data analysis for the paper:
Gloriana Chaverri, Maria Sagot, Jennifer L. Stynoski, Marcelo Araya-Salas, Yimen Araya-Ajoy, Martina Nagy, Mirjam Knörnschild, Silvia Chaves-Ramírez, Nicole Rose, Mariela Sánchez-Chavarría, Yelananie Jiménez-Torres, Deilyn Ulloa-Sanabria, Hellen Solís-Hernández, Gerald G. Carter. In review. Contact-calling rates within groups of disc-winged bats vary by caller but not by the receiver’s identity, kinship, or association. Philosophical Transactions of the Royal Society B.
Source code and data found at https://github.com/morceglo/Vocal-interactions-Thyroptera-tricolor
Load packages
Code
# load function from sketchy
source("https://raw.githubusercontent.com/maRce10/sketchy/main/R/load_packages.R")
# install/ load packages
sketchy::load_packages(packages = c("igraph", "bipartite", "asnipe",
"assortnet", "ggplot2", "ggmap", "ecodist", "igraphdata", "statnet",
"RColorBrewer", "tidyverse", "geosphere", "mapproj", "sna", "stringr",
"lme4", "glmmTMB", "broom.mixed", "boot", "patchwork", "performance",
"MASS"))1 Visualize social networks
1.1 Prepare data
Code
# Using
# https://dshizuka.github.io/networkanalysis/example_make_sparrow_network.html
# Input data
batdat21 = read.csv("./data/raw/associations_2021.csv", header = TRUE,
sep = ";", colClasses = c("numeric", "character", "Date", "character",
"character", "character", "character", "character", "character",
"character", "character", "character", "character"))
batdat22 = read.csv("./data/raw/associations_2022.csv", header = TRUE,
sep = ";", colClasses = c("numeric", "character", "Date", "character",
"character", "character", "character", "character", "character",
"character", "character"))
# Add attribute data
attrib_21 <- read.csv("./data/raw/supp_21.csv", header = TRUE, sep = ";",
colClasses = c("character", "character"))
attrib_22 <- read.csv("./data/raw/supp_22.csv", header = TRUE, sep = ";",
colClasses = c("character", "character"))
# Prepare association data 2021
batcols = grep("Bat", colnames(batdat21))
# batcols batdat21[,batcols] gather(batdat21[,batcols])
bat.ids = unique(gather(batdat21[, batcols])$value)
# bat.ids
bat.ids = bat.ids[!is.na(bat.ids)]
# bat.ids
bat.ids_df_21 <- data.frame(bat.ids)
# csv_path <- 'bats_21.csv' write.csv(bat.ids_df_21, csv_path,
# row.names = FALSE)
m1_21 = apply(batdat21[, batcols], 1, function(x) match(bat.ids, x))
# m1_21
m1_21[is.na(m1_21)] = 0
m1_21[m1_21 > 0] = 1
# m1_21
rownames(m1_21) = bat.ids #rows are bat ids
colnames(m1_21) = paste("group", 1:ncol(m1_21), sep = "_") #columns are group IDs (just 'group_#')
# m1_21
# rowSums(m1_21) set.seed(2) png(file='./output/Histogram
# captures 2021.png', width=600, height=400) par(oma=c(1,1,1,1))
# # all sides have 3 lines of space par(mar=c(5,4,4,2) + 0.8)
# hist(rowSums(m1_21), main='', xlab='Number of observations',
# ylab='Frequency', cex.lab = 2.5, cex.axis = 1.5) dev.off()
average_n_21 <- mean(rowSums(m1_21)) #average number of times a single bat was observed
# average_n_21
sd_n_21 <- sd(rowSums(m1_21)) #sd number of times a single bat was observed
# sd_n_21
average_gs_21 <- mean(colSums(m1_21)) #average group size
# average_gs_21
sd_gs_21 <- sd(colSums(m1_21)) #sd group size
# sd_gs_21
# Prepare association data 2022
batcols = grep("Bat", colnames(batdat22))
# batcols batdat22[,batcols] gather(batdat22[,batcols])
bat.ids = unique(gather(batdat22[, batcols])$value)
# bat.ids
bat.ids = bat.ids[is.na(bat.ids) == F]
# bat.ids
bat.ids_df_22 <- data.frame(bat.ids)
# csv_path = 'bats_22.csv' write.csv(bat.ids_df_22, csv_path,
# row.names = FALSE)
m1_22 = apply(batdat22[, batcols], 1, function(x) match(bat.ids, x))
# m1_22
m1_22[is.na(m1_22)] = 0
m1_22[m1_22 > 0] = 1
# m1_22
rownames(m1_22) = bat.ids #rows are bat ids
colnames(m1_22) = paste("group", 1:ncol(m1_22), sep = "_") #columns are group IDs (just 'group_#')
# m1_22
# rowSums(m1_22) set.seed(2) png(file='./output/Histogram
# captures 2022.png', width=600, height=400) par(oma=c(1,1,1,1))
# # all sides have 3 lines of space par(mar=c(5,4,4,2) + 0.8)
# hist (rowSums(m1_22), main='', xlab='Number of observations',
# ylab='Frequency', cex.lab = 2.5, cex.axis = 1.5) dev.off()
average_n_22 <- mean(rowSums(m1_22)) #average number of times a single bat was observed
# average_n_22
sd_n_22 <- sd(rowSums(m1_22)) #sd number of times a single bat was observed
# sd_n_22
average_gs_22 <- mean(colSums(m1_22)) #average group size
# average_gs_22
sd_gs_22 <- sd(colSums(m1_22)) #average group size
# sd_gs_221.2 Plot networks
Code
Generating 76 x 76 matrixCode
g_21=graph_from_adjacency_matrix(adj_21, "undirected", weighted=T) #the igraph object
V(g_21)$sex <- factor(attrib_21[match(V(g_21)$name, attrib_21$bat_id), "sex"])
mismatched_names <- setdiff(V(g_21)$name, attrib_21$bat_id) #check for issues
#Create adjacency data for associations (2022)
adj_22=get_network(t(m1_22), data_format="GBI", association_index = "SRI") # the adjacency matrixGenerating 76 x 76 matrixCode
g_22=graph_from_adjacency_matrix(adj_22, "undirected", weighted=T) #the igraph object
V(g_22)$sex <- factor(attrib_22[match(V(g_22)$name, attrib_22$bat_id), "sex"])
mismatched_names <- setdiff(V(g_22)$name, attrib_22$bat_id) #check for issues
#Select individuals for which more than 10 observations were made (not incorporated into adjacency data for associations, above)
#m2_21=m1_21[which(rowSums(m1_21)>10),]
#adj_21=get_network(t(m2_21), data_format="GBI","SRI")
#write.csv(adj_21, "assoc_21.csv", row.names=TRUE)
#rowSums(m2_21)
#m2_22=m1_22[which(rowSums(m1_22)>10),]
#adj_22=get_network(t(m2_22), data_format="GBI","SRI")
#rowSums(m2_22)
#write.csv(adj_22, "assoc_22.csv", row.names=TRUE)
#rowSums(m2_22)
#Create communities (2021)
com_21=fastgreedy.community(g_21) #community detection method
# com_21[[11]]
# Assign community names
mapping <- c(
11, # Old membership 1 should be corrected to 11
8, # Old membership 2 should be corrected to 8
1, # Old membership 3 should be corrected to 1
2, # Old membership 4 should be corrected to 2
6, # Old membership 5 should be corrected to 6
9, # Old membership 6 should be corrected to 9
4, # Old membership 7 should be corrected to 4
10, # Old membership 8 should be corrected to 10
7, # Old membership 9 should be corrected to 7
5, # Old membership 10 should be corrected to 5
3 # Old membership 11 should be corrected to 3
)
corrected_memberships <- mapping[com_21$membership]
com_21$membership <- corrected_memberships
community_assignments <- com_21$membership
color_mapping <- c(
"female" = "slateblue",
"male" = "gold",
"unknown" = "gray"
)
node_sex <- V(g_21)$sex
node_colors <- sapply(node_sex, function(sex) {
color_mapping[sex]
})
V(g_21)$color <- node_colors
#Prepare to save figure
# png(file="./output/Association networks 2021.png",
# width=600, height=600)
par(oma=c(1,1,1,1), mar=c(5,4, 4, 2) + 0.8, mfrow = c(1, 2))
# Iterate through com_21 and assign community names to nodes in g_21
plot(
g_21,
vertex.color = node_colors, # Use the calculated node colors
vertex.size = 10, # Adjust the size of the nodes as needed
vertex.label = community_assignments, # Use the extracted community assignments as labels
edge.width=E(g_21)$weight*10
)
title(main = "2021", cex.main = 2)
legend("topleft", legend = names(color_mapping), fill = color_mapping, cex = 1.5)
#Create communities (2022)
com_22=fastgreedy.community(g_22) #community detection method
# Assign community names
mapping <- c(
3,
4,
1,
9,
10,
11,
7,
5,
6,
2,
12,
8
)
corrected_memberships <- mapping[com_22$membership]
com_22$membership <- corrected_memberships
community_assignments <- com_22$membership
color_mapping <- c(
"female" = "slateblue",
"male" = "gold",
"unknown" = "gray"
)
node_sex <- V(g_22)$sex
node_colors <- sapply(node_sex, function(sex) {
color_mapping[sex]
})
V(g_22)$color <- node_colors
#Prepare to save figure
# png(file="./output/Association networks 2022.png",width=600, height=600)
# Iterate through com_21 and assign community names to nodes in g_21
plot(
g_22,
vertex.color = node_colors, # Use the calculated node colors
vertex.size = 10, # Adjust the size of the nodes as needed
vertex.label = community_assignments, # Use the extracted community assignments as labels
edge.width=E(g_22)$weight*10
)
title(main = "2022", cex.main = 2)
1.3 Estimate modulatity
Modularity:
- 2021: 0.87
- 2022: 0.9
2 Calling rate repeatability
2.1 Negative binomial model
Code
# Data manipulation Joining the two data sets
d1 <- read.csv("./data/raw/vocal_interactions_2021.csv", sep = ";")
d2 <- read.csv("./data/raw/vocal_interactions_2022.csv", sep = ";")
colnames(d2)[6] <- "Dyad"
d <- rbind(d1, d2)
## Adding observation level random effect
d$ob <- 1:nrow(d)
# Stats models Model for inquiry
mod_inquiry <- glmer.nb(n_inquiry ~ 1 + (1 | Bat_roosting) + (1 |
Bat_flying) + (1 | Group), data = d)
summary(mod_inquiry)Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: Negative Binomial(2.7219) ( log )
Formula: n_inquiry ~ 1 + (1 | Bat_roosting) + (1 | Bat_flying) + (1 |
Group)
Data: d
AIC BIC logLik deviance df.resid
4596.7 4617.2 -2293.3 4586.7 446
Scaled residuals:
Min 1Q Median 3Q Max
-1.62529 -0.48778 0.07217 0.47606 2.82516
Random effects:
Groups Name Variance Std.Dev.
Bat_roosting (Intercept) 8.029e-11 8.961e-06
Bat_flying (Intercept) 2.348e-01 4.845e-01
Group (Intercept) 2.141e-02 1.463e-01
Number of obs: 451, groups: Bat_roosting, 74; Bat_flying, 73; Group, 14
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.12062 0.07845 52.53 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
optimizer (Nelder_Mead) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')Code
### Extracting variance components
mu_l <- fixef(mod_inquiry)[1]
mu <- exp(mu_l)
theta <- getME(mod_inquiry, "glmer.nb.theta")
VBF_I_l <- VarCorr(mod_inquiry)$Bat_flying[1, 1] #Variance Bats Flying
VBF_I <- (exp(VBF_I_l) - 1) * exp(2 * mu_l + VBF_I_l)
VBR_I_l <- VarCorr(mod_inquiry)$Bat_roosting[1, 1] #Variance Bats Roosting
VBR_I <- (exp(VBR_I_l) - 1) * exp(2 * mu_l + VBR_I_l)
VBG_I_l <- VarCorr(mod_inquiry)$Group[1, 1] #Variance GRoup
VBG_I <- (exp(VBG_I_l) - 1) * exp(2 * mu_l + VBG_I_l)
NBV_I <- mu + mu^2/theta
#### Repeatabilty roosting individual
VBR_I/(VBR_I + VBF_I + NBV_I + VBG_I) (Intercept)
1.084486e-10 (Intercept)
0.465914 Code
Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: Negative Binomial(4.214) ( log )
Formula: n_inquiry ~ 1 + (1 | Bat_roosting) + (1 | Bat_flying) + (1 |
Group)
Data: d2
AIC BIC logLik deviance df.resid
2606.7 2624.4 -1298.4 2596.7 250
Scaled residuals:
Min 1Q Median 3Q Max
-1.93767 -0.52024 0.04805 0.49889 2.73912
Random effects:
Groups Name Variance Std.Dev.
Bat_roosting (Intercept) 7.930e-13 8.905e-07
Bat_flying (Intercept) 2.221e-01 4.712e-01
Group (Intercept) 3.879e-02 1.970e-01
Number of obs: 255, groups: Bat_roosting, 71; Bat_flying, 69; Group, 14
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.21595 0.08999 46.85 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
optimizer (Nelder_Mead) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')Code
### Extracting variance components
mu_l <- fixef(mod_inquiry)[1]
mu <- exp(mu_l)
theta <- getME(mod_inquiry, "glmer.nb.theta")
VBF_I_l <- VarCorr(mod_inquiry)$Bat_flying[1, 1] #Variance Bats Flying
VBF_I <- (exp(VBF_I_l) - 1) * exp(2 * mu_l + VBF_I_l)
VBR_I_l <- VarCorr(mod_inquiry)$Bat_roosting[1, 1] #Variance Bats Roosting
VBR_I <- (exp(VBR_I_l) - 1) * exp(2 * mu_l + VBR_I_l)
VBG_I_l <- VarCorr(mod_inquiry)$Group[1, 1] #Variance GRoup
VBG_I <- (exp(VBG_I_l) - 1) * exp(2 * mu_l + VBG_I_l)
NBV_I <- mu + mu^2/theta
#### Repeatabilty roosting individual
VBR_I/(VBR_I + VBF_I + NBV_I) (Intercept)
1.409571e-12 (Intercept)
0.5519081 Code
Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: Negative Binomial(0.1622) ( log )
Formula: n_response ~ 1 + (1 | Bat_roosting) + (1 | Bat_flying) + (1 |
Group)
Data: d
AIC BIC logLik deviance df.resid
3742.3 3762.8 -1866.1 3732.3 446
Scaled residuals:
Min 1Q Median 3Q Max
-0.4026 -0.4017 -0.3616 0.2471 4.1603
Random effects:
Groups Name Variance Std.Dev.
Bat_roosting (Intercept) 1.118e+00 1.057e+00
Bat_flying (Intercept) 5.331e-13 7.301e-07
Group (Intercept) 1.489e-12 1.220e-06
Number of obs: 451, groups: Bat_roosting, 74; Bat_flying, 73; Group, 14
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.074 0.213 19.12 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
optimizer (Nelder_Mead) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')Code
### Extracting variance components
mu_l <- fixef(mod_response)[1] #+ fixef(mod_response)[2]*mean(d$n_inquiry, na.rm=TRUE)
mu <- exp(mu_l)
theta <- getME(mod_response, "glmer.nb.theta")
VBF_R_l <- VarCorr(mod_response)$Bat_flying[1, 1] #Variance Bats Flying
VBF_R <- (exp(VBF_R_l) - 1) * exp(2 * mu_l + VBF_R_l)
VBR_R_l <- VarCorr(mod_response)$Bat_roosting[1, 1] #Variance Bats Roosting
VBR_R <- (exp(VBR_R_l) - 1) * exp(2 * mu_l + VBR_R_l)
VBG_R_l <- VarCorr(mod_response)$Group[1, 1] #Variance GRoup
VBG_R <- (exp(VBG_R_l) - 1) * exp(2 * mu_l + VBG_R_l)
NBV_R <- mu + mu^2/theta
#### Repeatabilty roosting individual
VBR_R/(VBR_R + VBF_R + NBV_R + VBG_R)(Intercept)
0.5044 (Intercept)
4.274015e-14 2.2 Poisson model
Code
# #Data manipulation ## Joining the two data sets
# d1<-read.csv('vocal_interactions_2021.csv', sep=';')
# d2<-read.csv('vocal_interactions_2022.csv', sep=';')
# colnames(d2)[6]<-'Dyad' d<-rbind(d1,d2)
## Adding observation level random effect
d$ob <- 1:nrow(d)
# Stats models Model for inquiry with all observations
mod_inquiry <- glmer(n_inquiry ~ 1 + (1 | Bat_roosting) + (1 | Bat_flying) +
(1 | Group) + (1 | ob), data = d, family = "poisson")
### Extracting variance components
VBF_I <- VarCorr(mod_inquiry)$Bat_flying[1, 1] #Variance Bats Flying
VBR_I <- VarCorr(mod_inquiry)$Bat_roosting[1, 1] #Variance Bats Roosting
VBG_I <- VarCorr(mod_inquiry)$Group[1, 1] #Variance Group
VOD_I <- VarCorr(mod_inquiry)$ob[1, 1] #Variance Overddispersion
PV_I <- log(1/exp(fixef(mod_inquiry)[1]) + 1) #Variance poisson processs
#### Repeatabilty roosting individual
VBR_I/(VBR_I + VOD_I + VBF_I + PV_I + VBG_I)(Intercept)
0.07582562 (Intercept)
0.3913205 Code
## Model for inquiry only when someone responded
d2 <- d[d$n_response > 0, ]
mod_inquiry2 <- glmer(n_inquiry ~ 1 + (1 | Bat_roosting) + (1 | Bat_flying) +
(1 | ob), data = d2, family = "poisson")
### Extracting variance components
VBF_I <- VarCorr(mod_inquiry2)$Bat_flying[1, 1] #Variance Bats Flying
VBR_I <- VarCorr(mod_inquiry2)$Bat_roosting[1, 1] #Variance Bats Roosting
VBG_I <- VarCorr(mod_inquiry)$Group[1, 1]
VOD_I <- VarCorr(mod_inquiry2)$ob[1, 1] #Variance Overddispersion
PV_I <- log(1/exp(fixef(mod_inquiry2)[1]) + 1) #Variance poisson processs
#### Repeatabilty roosting individual
VBR_I/(VBR_I + VOD_I + VBF_I + PV_I + VBG_I)(Intercept)
0.02516122 (Intercept)
0.5240666 Code
Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: poisson ( log )
Formula: n_response ~ 1 + (1 | Bat_roosting) + (1 | Bat_flying) + (1 |
Group) + (1 | ob)
Data: d
AIC BIC logLik deviance df.resid
3816.4 3837.0 -1903.2 3806.4 446
Scaled residuals:
Min 1Q Median 3Q Max
-0.75743 -0.39416 -0.00652 0.01611 0.17248
Random effects:
Groups Name Variance Std.Dev.
ob (Intercept) 8.3383 2.8876
Bat_roosting (Intercept) 7.8244 2.7972
Bat_flying (Intercept) 0.4551 0.6746
Group (Intercept) 1.0819 1.0401
Number of obs: 451, groups:
ob, 451; Bat_roosting, 74; Bat_flying, 73; Group, 14
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.1595061 0.0004604 2518 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
optimizer (Nelder_Mead) convergence code: 0 (OK)
Model failed to converge with max|grad| = 0.15353 (tol = 0.002, component 1)
Model is nearly unidentifiable: very large eigenvalue
- Rescale variables?Code
### Extracting variance components
VBF_R <- VarCorr(mod_response)$Bat_flying[1, 1] #Variance Bats Flying
VBR_R <- VarCorr(mod_response)$Bat_roosting[1, 1] #Variance Bats Roosting
VBG_R <- VarCorr(mod_response)$Group[1, 1]
VOD_R <- VarCorr(mod_response)$ob[1, 1] #Variance Overddispersion
PV_R <- log(1/exp(fixef(mod_response)[1]) + 1) #Variance poisson processs
## Repeatabilty roosting individual
VBR_R/(VBR_R + VOD_R + VBF_R + PV_R + VBG_R)(Intercept)
0.4353549 (Intercept)
0.02532077 3 Calling, association and kinship
3.1 Prepare data
Code
# Code by Gerald Carter, gcarter1640@gmail.com
### clean and wrangle data
# create function to convert matrix to data-frame
matrix_to_df <- function(m1){
data.frame(dyad = paste(rownames(m1)[col(m1)], colnames(m1)[row(m1)], sep="_"),
value = c(t(m1)), stringsAsFactors = FALSE)
}
# get kinship data
k <- as.matrix(read.csv('./data/raw/KinshipMatrix.csv', sep= ",", row.names = 1))
# check symmetry
mean(k==t(k), na.rm=T)
# get sex data
sex1 <- read.csv('supp_21.csv', sep= ";")
sex2 <- read.csv('supp_22.csv', sep= ";")
# get association data
a21 <- read.csv('associations_2021.csv', sep= ";")
a22 <- read.csv('associations_2022.csv', sep= ";")
# get calling data
i21 <- read.csv('vocal_interactions_2021.csv', sep= ";")
i22 <- read.csv('vocal_interactions_2022.csv', sep= ";")
# tidy sex data
sex <-
rbind(sex1,sex2) %>%
mutate(bat= paste0('bat', bat_id)) %>%
group_by(bat, sex) %>%
summarize(n= n()) %>%
dplyr::select(bat, sex) %>%
ungroup()
# tidy 2021 association data
assoc21 <-
a21 %>%
pivot_longer(Bat1:Bat10, names_to = 'names', values_to = "bat") %>%
filter(!is.na(bat)) %>%
mutate(bat= paste0('bat', bat)) %>%
mutate(group= paste0('group', Group)) %>%
mutate(date= dmy(Date)) %>%
dplyr::select(obs= data_id, group, date, bat)
# tidy 2022 association data
assoc22 <-
a22 %>%
pivot_longer(Bat1:Bat8, names_to = 'names', values_to = "bat") %>%
filter(!is.na(bat)) %>%
mutate(bat= paste0('bat', bat)) %>%
mutate(group= paste0('group', Group)) %>%
mutate(date= dmy(Date)) %>%
dplyr::select(obs= data_id, group, date, bat)
# combine association data
assoc <-
rbind(assoc21, assoc22) %>%
mutate(obs= paste(obs,group,date, sep= "_")) %>%
dplyr::select(bat, obs) %>%
as.data.frame()
# get number of observations per bat
obs_per_bat <-
assoc %>%
group_by(bat) %>%
summarize(n.obs=n()) %>%
arrange(desc(n.obs))
# range is 1 to 51 times
range(obs_per_bat$n.obs)
# pick threshold of observations for including bats in analysis
# if they have fewer observations than we make their association rates NA
threshold <- 4
# plot number of observations per bat
assoc %>%
group_by(bat) %>%
summarize(n.obs=n()) %>%
ggplot(aes(x=n.obs))+
geom_histogram(fill="light blue", color="black")+
geom_vline(xintercept= threshold, color= 'red', size=1)+
xlab("number of observations of bat")+
ylab("count of bats")
# get bats seen fewer than threshold number
low.n.bats <-
assoc %>%
group_by(bat) %>%
summarize(n.obs=n()) %>%
filter(n.obs < threshold) %>%
pull(bat)
low.n.bats
length(low.n.bats)
# 14 bats were seen fewer than 4 times
### get SRI from asnipe
# get association rates as group-by-individual
gbi <- get_group_by_individual(assoc, data_format="individuals")
# get association network of SRI
assoc.net<- get_network(association_data=gbi, data_format = "GBI", association_index = "SRI")
# check symmetry
mean(assoc.net==t(assoc.net))
# get SRI values for every undirected pair as dataframe
SRI <- matrix_to_df(assoc.net)
# tidy kinship
kinship <-
k %>%
matrix_to_df() %>%
separate(dyad, into= c("bat1", "bat2")) %>%
mutate(bat2= str_remove(bat2, "X")) %>%
mutate(bat1= paste0('bat', bat1), bat2= paste0('bat', bat2)) %>%
filter(bat1!=bat2) %>%
# label undirected pair
mutate(dyad= ifelse(bat1<bat2, paste(bat1,bat2, sep="_"), paste(bat2,bat1, sep="_"))) %>%
group_by(dyad) %>%
summarize(kinship= first(value))
# count trials
nrow(i21)
nrow(i22)
# trials with missing data
rbind(i21,i22) %>%
as_tibble() %>%
filter(is.na(n_response)) %>%
nrow()
#2
# trials where neither bat called
rbind(i21,i22) %>%
as_tibble() %>%
filter(n_inquiry==0) %>%
nrow()
# 5
# fix and tidy calling data
calling <-
rbind(i21,i22) %>%
mutate(Date= as.character(mdy(Date))) %>%
mutate(group= paste(substr(Date, start=1, stop=4),Group, sep="_")) %>%
# need to recreate these labels because messed up by excel in raw data
separate(Dyad, into=c("bat_flying", "bat_roosting"), remove=F) %>%
mutate(bat_flying= paste0('bat', bat_flying), bat_roosting= paste0('bat', bat_roosting)) %>%
# label undirected dyads
mutate(dyad= ifelse(bat_flying<bat_roosting, paste(bat_flying, bat_roosting, sep= "_"), paste(bat_roosting, bat_flying, sep= "_"))) %>%
dplyr::select(date= Date, group, bat_flying, bat_roosting, dyad, flight_time, n_inquiry, n_response)
# combine all data
d <-
# start with calling data
calling %>%
# add SRI
mutate(sri= SRI$value[match(.$dyad, SRI$dyad)]) %>%
# add kinship
mutate(kinship= kinship$kinship[match(.$dyad, kinship$dyad)]) %>%
# add sexes of flying and roosting bats
mutate(flying.sex= sex$sex[match(.$bat_flying, sex$bat)]) %>%
mutate(roosting.sex= sex$sex[match(.$bat_roosting, sex$bat)]) %>%
# label sex combination for flying-->roosting dyad
mutate(dyad.sex= paste(flying.sex, roosting.sex, sep= "-->")) %>%
# label sex combination for dyad (male, female, both)
mutate(udyad.sex = case_when(
flying.sex == 'female' & roosting.sex == 'female' ~ "female-female",
flying.sex == 'male' & roosting.sex == 'male' ~ "male-male",
TRUE ~ "mixed-sex")) %>%
# replace zero sri (never previously seen together) with NA because association is not 0
mutate(sri = if_else(sri==0, NA, sri)) %>%
# if flying bat seen fewer than 4 times, then SRI is NA
mutate(sri = if_else(bat_flying %in% low.n.bats, NA, sri)) %>%
# if roosting bat seen fewer than 4 times, then SRI is NA
mutate(sri = if_else(bat_roosting %in% low.n.bats, NA, sri))
# inspect and fix cases where flying bats in more than 2 groups
t <-
d %>%
group_by(bat_flying, date, group) %>%
summarize(n=n()) %>%
arrange(date) %>%
group_by(bat_flying, group) %>%
summarize(n=n()) %>%
group_by(bat_flying) %>%
summarize(n=n()) %>%
filter(n>2) %>%
pull(bat_flying)
d %>%
dplyr::select(date, group, bat_flying) %>%
filter(bat_flying %in% t) %>%
arrange(bat_flying) %>%
as.data.frame()
# some of these seem impossible and must be errors
# fix impossible group labels
d$group[which(d$date=="2021-07-03" & d$bat_flying=="bat982126051278521")] <- "2021_9"
d$group[which(d$date=="2021-07-03" & d$bat_flying=="bat982126058484300")] <- "2021_9"
# fix cases where roosting bats in more than 2 groups
t <-
d %>%
group_by(bat_roosting, date, group) %>%
summarize(n=n()) %>%
arrange(date) %>%
group_by(bat_roosting, group) %>%
summarize(n=n()) %>%
group_by(bat_roosting) %>%
summarize(n=n()) %>%
filter(n>2) %>%
pull(bat_roosting)
d %>%
dplyr::select(date, group, bat_roosting) %>%
filter(bat_roosting %in% t) %>%
arrange(bat_roosting) %>%
as.data.frame()
# fix group labels
d$group[which(d$date=="2021-07-03" & d$bat_roosting=="bat982126058484300")] <- "2021_9"
# group 9 split into groups 9,1 and 9,2 during 2022
# There are interesting movements between group 9 and 10 and between groups 9,1 and 9,2
# For this analysis, I'm going to consider group 9,1 and 9,2 as the same group
d$group[which(d$group=="2022_9,1" | d$group=="2022_9,2")] <- "2021_9"
# remove trials without inquiry or response calls
d <-
d %>% # 453 rows
as_tibble() %>%
filter(! (is.na(n_inquiry) & is.na(n_response))) %>% # 451 rows
filter(n_inquiry>0) %>% # 446 rows
mutate(year= substr(date, 1,4)) # add year
# save data as csv
write.csv(d, file = paste0('./processed/calling-association-kinship-data.csv'), row.names= F)3.2 Explore data
Code
# set ggplot theme
theme_set(theme_bw(base_size = 14))
# get data
raw <- read.csv("./data/processed/calling-association-kinship-data.csv")
###### Create functions to estimate confindence intervals
###### function to get mean and 95% CI of values x via
###### bootstrapping
boot_ci <- function(x, perms = 5000, bca = F) {
get_mean <- function(x, d) {
return(mean(x[d]))
}
x <- as.vector(na.omit(x))
mean <- mean(x)
if (bca) {
boot <- boot.ci(boot(data = x, statistic = get_mean, R = perms,
parallel = "multicore", ncpus = 4), type = "bca")
low <- boot$bca[1, 4]
high <- boot$bca[1, 5]
} else {
boot <- boot.ci(boot(data = x, statistic = get_mean, R = perms,
parallel = "multicore", ncpus = 4), type = "perc")
low <- boot$perc[1, 4]
high <- boot$perc[1, 5]
}
c(low = low, mean = mean, high = high, N = round(length(x)))
}
# function to get mean and 95% CI via bootstrapping of values y
# within grouping variable x
boot_ci2 <- function(d = d, y = d$y, x = d$x, perms = 5000, bca = F) {
df <- data.frame(effect = unique(x))
df$low <- NA
df$mean <- NA
df$high <- NA
df$n.obs <- NA
for (i in 1:nrow(df)) {
ys <- y[which(x == df$effect[i])]
if (length(ys) > 1 & var(ys) > 0) {
b <- boot_ci(y[which(x == df$effect[i])], perms = perms,
bca = bca)
df$low[i] <- b[1]
df$mean[i] <- b[2]
df$high[i] <- b[3]
df$n.obs[i] <- b[4]
} else {
df$low[i] <- min(ys)
df$mean[i] <- mean(ys)
df$high[i] <- max(ys)
df$n.obs[i] <- length(ys)
}
}
df
}
# look at data head(raw) group is the year and social group dyad
# is the undirected pair (flying and roosting bat) in
# alphanumeric order sri is simple ratio index of association
# flying.sex is sex of flying bat roosting.sex is sex of
# roosting bat dyad.sex is the sex of the flying and roosting
# bat udyad.sex is females, males, or mixed
# add a few variables to the data
d <- raw %>%
# get log count of inquiry calls
mutate(log_inquiry = log(n_inquiry)) %>%
# rescale kinship so 1 unit is 0.5
mutate(kinship2 = kinship * 2) %>%
# if the roosting bat leaves the roost (escapes) then flight
# time is < 300 s
mutate(escape = as.numeric(flight_time < 300)) %>%
# convert flight time from seconds to minutes (for easier
# interpretation)
mutate(flight_time2 = flight_time/60)
###### Describe the data-------------
# how many trials? d %>% nrow() 446
# # how many undirected pairs have vocal data? d %>% pull(dyad)
# %>% unique() %>% length 139 undirected pairs
# d %>% group_by(bat_flying, bat_roosting) %>% summarize(n=n())
# 254 directed pairs
# how many groups? n_distinct(d$group) 23
# how many of these have association data d %>% group_by(dyad)
# %>% summarize(sri= mean(sri)) %>% filter(!is.na(sri)) 133
# pairs have association
# how many of these have kinship data d %>% group_by(dyad) %>%
# summarize(kinship= mean(kinship)) %>% filter(!is.na(kinship))
# 82 pairs have kinship data
# how many related undirected pairs d %>% group_by(dyad) %>%
# summarize(kinship= mean(kinship)) %>% filter(kinship>0) 32 kin
# pairs
# how many unrelated undirected pairs d %>% group_by(dyad) %>%
# summarize(kinship= mean(kinship)) %>% filter(kinship==0) 50
# nonkin pairs
# what is mean kinship in group?
grp.kin <- d %>%
group_by(dyad, group) %>%
summarize(kinship = mean(kinship, na.rm = T)) %>%
group_by(group) %>%
summarize(kinship = mean(kinship, na.rm = T), n = n()) %>%
as.data.frame()
set.seed(123)
# grp.kin %>% pull(kinship) %>% boot_ci(bca=T) 0.23, 95% CI =
# [0.16, 0.31]
# how many groups have kin? grp.kin %>% filter(kinship>0) 17
# how many groups have no kin? grp.kin %>% filter(kinship==0) 4
# plot distribution of kinship
d %>%
ggplot(aes(x = kinship)) + facet_wrap(~udyad.sex, ncol = 1, scales = "free_y") +
geom_histogram(fill = "light blue", color = "black") + ggtitle("pairwise kinship")
Code
# # check categories d %>% filter(!is.na(kinship)) %>%
# group_by(dyad) %>% summarize(kinship= mean(kinship)) %>%
# group_by(kinship) %>% summarize(n=n())
##### Effect of kinship on association-------------
# plot distribution of assoc
d %>%
ggplot(aes(x = sri)) + facet_wrap(~udyad.sex, ncol = 1, scale = "free_y") +
geom_histogram(fill = "light blue", color = "black") + xlab("association rate") +
ggtitle("pairwise SRI values")
Code
# looks ok
# what is mean and 95% CI of assoc among flying and roosting
# bats? set.seed(123) d %>% group_by(dyad) %>% summarize(sri =
# mean(sri)) %>% pull(sri) %>% boot_ci(bca=T) 0.51, 95% CI =
# [0.47, 0.55] a bit low, expected from past work is 0.76
# now only nonkin
set.seed(123)
k1 <- d %>%
filter(kinship == 0) %>%
group_by(dyad) %>%
summarize(sri = mean(sri)) %>%
pull(sri) %>%
boot_ci(bca = T) %>%
c(kinship = 0)
# k1 0.572, 95% CI = [0.500, 0.636]
# now only kin set.seed(123) d %>% filter(kinship>0) %>%
# group_by(dyad) %>% summarize(sri = mean(sri)) %>% pull(sri)
# %>% boot_ci(bca=T) 0.686, 95% CI = [0.631, 0.741]
# now only close kin
set.seed(123)
k2 <- d %>%
filter(kinship == 0.5) %>%
group_by(dyad) %>%
summarize(sri = mean(sri)) %>%
pull(sri) %>%
boot_ci(bca = T) %>%
c(kinship = 0.5)
# k2 0.687, 95% CI = [0.628, 0.750]
# now only 0.25 kin
set.seed(123)
k3 <- d %>%
filter(kinship == 0.25) %>%
group_by(dyad) %>%
summarize(sri = mean(sri)) %>%
pull(sri) %>%
boot_ci(bca = T) %>%
c(kinship = 0.25)
# k3 0.684, 95% CI = [0.542, 0.783]
# compile means
k <- rbind(k1, k2, k3) %>%
data.frame() %>%
mutate(kin = as.character(kinship)) %>%
mutate(kin2 = kinship > 0) %>%
mutate(assoc = mean)
# plot association by kinship
(kin.dyads.plot <- d %>%
filter(!is.na(kinship)) %>%
group_by(dyad) %>%
summarize(kinship = mean(kinship), assoc = mean(sri), udyad.sex = first(udyad.sex)) %>%
mutate(kin2 = kinship > 0) %>%
mutate(kin = as.character(kinship)) %>%
ggplot(aes(x = kin, y = assoc, group = kin, color = kin2)) + geom_jitter(size = 2,
width = 0.1, height = 0, aes(shape = udyad.sex)) + geom_point(data = k,
position = position_nudge(x = 0.25), size = 3) + geom_errorbar(data = k,
aes(ymin = low, ymax = high, width = 0.1), position = position_nudge(x = 0.25),
size = 1) + xlab("kinship") + ylab("association rate (simple ratio index)") +
scale_color_manual(values = c("darkgrey", "darkred")) + theme(legend.position = "none",
legend.title = element_blank()))
Code
# save as PDF ggsave( 'kin_association.pdf', plot =
# kin.dyads.plot, scale = 1, width = 3, height = 5, units =
# 'in', dpi = 300)
##### Effect of flight time on calling ----------
### How does flight time and count of inquiry calls predict
### count of response calls?
# plot number of response calls by flight time
t1 <- d %>%
mutate(kinship = kinship > 0.1) %>%
ggplot(aes(x = flight_time, y = n_response)) + geom_point(size = 2,
alpha = 0.3, aes(color = kinship, shape = kinship)) + geom_smooth(method = "glm.nb") +
scale_color_manual(values = c("darkgrey", "darkred")) + xlab("flight time (seconds)") +
ylab("response call count") + theme(legend.position = "none")
# plot number of inquiry calls by flight time
t2 <- d %>%
mutate(kinship = kinship > 0.1) %>%
ggplot(aes(x = flight_time, y = n_inquiry)) + geom_point(size = 2,
alpha = 0.3, aes(color = kinship, shape = kinship)) + geom_smooth(method = "glm.nb") +
scale_color_manual(values = c("darkgrey", "darkred")) + xlab("flight time (seconds)") +
ylab("inquiry call count") + theme(legend.position = "none")
# plot together
t <- t1 + t2
# ggsave( filename= 'flight_time.pdf', plot = t, scale = 1,
# width = 7, height = 7, units = c('in', 'cm', 'mm', 'px'), dpi
# = 300)3.3 Statistical analysis
Code
# Overdispersion test
dispersion ratio = 6.288
Pearson's Chi-Squared = 2773.032
p-value = < 0.001Code
# Overdispersion test
dispersion ratio = 0.786
Pearson's Chi-Squared = 345.749
p-value = 1Code
# A tibble: 2 × 9
effect component term estimate std.error statistic p.value conf.low
<chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 fixed cond (Intercept) 11.9 1.51 19.5 7.77e-85 9.30
2 fixed cond flight_time2 1.43 0.0336 15.1 2.60e-51 1.36
# ℹ 1 more variable: conf.high <dbl>Code
# for every minute of flight time, the inquiry call count
# increases by a factor of 1.43 [1.36, 1.49]
# try poisson model for effect of flight time on response
# calling
t <- glmer(n_response ~ flight_time2 + log_inquiry + (1 | bat_flying) +
(1 | bat_roosting) + (1 | group), data = d, family = poisson)
check_overdispersion(t)# Overdispersion test
dispersion ratio = 63.477
Pearson's Chi-Squared = 27930.076
p-value = < 0.001Code
# Overdispersion test
dispersion ratio = 0.582
Pearson's Chi-Squared = 255.097
p-value = 1# Check for zero-inflation
Observed zeros: 191
Predicted zeros: 30
Ratio: 0.16Code
# A tibble: 4 × 9
effect component term estimate std.error statistic p.value conf.low conf.high
<chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 fixed cond (Int… 3.00 1.60 2.06 3.95e-2 1.08 8.91
2 fixed cond flig… 1.58 0.173 4.18 2.98e-5 1.26 1.95
3 fixed cond log_… 1.45 0.188 2.83 4.62e-3 1.12 1.87
4 fixed zi (Int… 0.644 0.0722 -3.93 8.65e-5 1.08 8.91Code
# for every minute of flight time, the response call count
# increases by a factor of 1.58 [1.26, 1.95]
##### Effect of inquiry calling on response calling------------
# plot with log counts
d %>%
ggplot(aes(x = log_inquiry, y = log(n_response + 1))) + geom_point(size = 2) +
geom_smooth(method = "lm") + xlab("log inquiry call count") +
ylab("log (x+1) response call count")
Code
# plot negative binomial curve plot all bats
(t1 <- d %>%
mutate(label = "all pairs") %>%
mutate(escape = as.logical(escape)) %>%
ggplot(aes(x = log_inquiry, y = n_response)) + facet_wrap(~label) +
geom_point(size = 2, alpha = 0.5) + geom_smooth(method = "glm.nb") +
xlab("log of inquiry call count") + ylab("response call count")) +
theme(legend.position = "none")
Code
# plot with kinship
(t2 <- d %>%
filter(!is.na(kinship)) %>%
mutate(escape = as.logical(escape)) %>%
mutate(kin = if_else(kinship > 0.1, "kin", "nonkin")) %>%
mutate(kin = factor(kin, levels = c("nonkin", "kin"))) %>%
ggplot(aes(x = log_inquiry, y = (n_response), color = kin, group = kin)) +
facet_wrap(~kin) + geom_point(size = 2, alpha = 0.5, aes(shape = escape)) +
geom_smooth(method = "glm.nb") + xlab("log of inquiry call count") +
ylab("response call count") + scale_color_manual(values = c("darkgrey",
"darkred")) + theme(legend.position = "none") + theme(axis.text.y = element_blank(),
axis.title.y = element_blank(), legend.position = "none"))
Code
# plot as linear model
(t3 <- d %>%
filter(!is.na(kinship)) %>%
mutate(kin = if_else(kinship > 0.1, "kin", "nonkin")) %>%
mutate(kin = factor(kin, levels = c("nonkin", "kin"))) %>%
ggplot(aes(x = log_inquiry, y = log(n_response + 1), color = kin,
group = kin)) + facet_wrap(~kin) + geom_point(size = 2, alpha = 0.5) +
geom_smooth(method = "lm") + xlab("log of inquiry call count") +
ylab("log of response call count") + scale_color_manual(values = c("darkgrey",
"darkred")) + theme(legend.position = "none"))
Code
# save as PDF ggsave( 'inquiry_response.pdf', plot =
# inquiry.response, width = 14, height = 7, units = 'in', dpi =
# 300)
# for plots, get residuals from negative binomial mixed effect
# model predicting number of responses by number of inquiry
# calls
fit <- glmmTMB(n_response ~ log_inquiry + offset(log(flight_time)) +
(1 | bat_flying) + (1 | bat_roosting) + (1 | group), data = d,
ziformula = ~1, family = nbinom2)
# get difference between observed response count and predicted
# response count add this response score to data for plotting
# this value represents response calling more than expected from
# inquiry calls
d$resid_response <- d$n_response - predict(fit, type = "response",
newdata = d, allow.new.levels = T)
# plot number of inquiry calls by kinship with partner
d %>%
filter(!is.na(kinship)) %>%
mutate(kin = as.character(kinship)) %>%
ggplot(aes(x = kin, y = n_inquiry, group = kin)) + geom_violin(fill = NA,
width = 1) + geom_jitter(size = 2, width = 0.1, height = 0, aes(color = udyad.sex)) +
xlab("kinship") + ylab("number of inquiry calls") + theme(legend.position = "top",
legend.title = element_blank())
Code
# plot number of response calls by kinship with partner
d %>%
filter(!is.na(kinship)) %>%
mutate(kin = as.character(kinship)) %>%
ggplot(aes(x = kin, y = n_response, group = kin)) + geom_violin(fill = NA,
width = 1) + geom_jitter(size = 2, width = 0.1, height = 0, aes(color = udyad.sex)) +
xlab("kinship") + ylab("number of response calls") + theme(legend.position = "top",
legend.title = element_blank())
Code
# plot residual response call variation
d %>%
filter(!is.na(kinship)) %>%
mutate(kin = as.character(kinship)) %>%
ggplot(aes(x = kin, y = resid_response, group = kin)) + geom_violin(fill = NA,
width = 1) + geom_jitter(size = 2, width = 0.1, height = 0, aes(color = udyad.sex)) +
xlab("kinship") + ylab("residual response call variation") + theme(legend.position = "top",
legend.title = element_blank())
Code
# create function to plot variable by kinship with means and 95%
# CIs
plot_kinship_effect <- function(d = d, y = d$sri, label = "label") {
# plot association by kinship
set.seed(123)
means <- d %>%
mutate(y = y) %>%
filter(!is.na(kinship)) %>%
group_by(dyad) %>%
summarize(kinship = mean(kinship), y = mean(y), udyad.sex = first(udyad.sex)) %>%
mutate(kin = ifelse(kinship > 0, "kin", "nonkin")) %>%
dplyr::select(kin, y) %>%
boot_ci2(y = .$y, x = .$kin, bca = T) %>%
rename(kin = effect) %>%
mutate(kin = factor(kin, levels = c("nonkin", "kin")))
points <- d %>%
mutate(y = y) %>%
filter(!is.na(kinship)) %>%
group_by(dyad) %>%
summarize(kinship = mean(kinship), y = mean(y), udyad.sex = first(udyad.sex)) %>%
mutate(kin = ifelse(kinship > 0, "kin", "nonkin")) %>%
mutate(kin = factor(kin, levels = c("nonkin", "kin")))
# plot means and 95% CI
means %>%
ggplot(aes(x = kin, y = mean, color = kin)) + geom_jitter(data = points,
aes(y = y), size = 1, alpha = 0.5, width = 0.1) + geom_point(position = position_nudge(x = 0.25),
size = 3) + geom_errorbar(aes(ymin = low, ymax = high, width = 0.1),
position = position_nudge(x = 0.25), size = 1) + xlab("") +
ylab(label) + scale_color_manual(values = c("darkgrey", "darkred")) +
coord_flip() + theme_bw() + theme(legend.position = "none")
}
# plot effect of kinship on SRI
(p1 <- plot_kinship_effect(d, d$sri, label = "association (SRI)"))
Code
Code
Code
Code
# plot number of inquiry calls by association and kinship (by
# dyad type)
t1 <- d %>%
filter(!is.na(kinship)) %>%
mutate(kin = ifelse(kinship > 0.1, "kin", "nonkin")) %>%
ggplot(aes(x = sri, y = log_inquiry)) + facet_wrap(~dyad.sex,
scales = "free_y") + geom_point(aes(color = kin), size = 2, alpha = 0.6) +
geom_smooth(method = "lm") + xlab("association rate (simple ratio index)") +
ylab("inquiry call count (log transformed)") + scale_colour_manual(values = c("darkred",
"grey")) + theme_bw() + theme(legend.position = "none", legend.title = element_blank())
# plot number of response calls by association and kinship (by
# dyad type)
t2 <- d %>%
filter(!is.na(kinship)) %>%
mutate(kin = ifelse(kinship > 0.1, "kin", "nonkin")) %>%
ggplot(aes(x = sri, y = log(n_response + 1))) + facet_wrap(~dyad.sex,
scales = "free_y") + geom_point(aes(color = kin), size = 2, alpha = 0.6) +
geom_smooth(method = "lm") + xlab("association rate (simple ratio index)") +
ylab("response call count (log transformed)") + scale_colour_manual(values = c("darkred",
"grey")) + theme_bw() + theme(legend.position = "none", legend.title = element_blank())
# plot residual response calling by association (by dyad type)-
# controlling for inquiry calls
t3 <- d %>%
filter(!is.na(kinship)) %>%
mutate(kin = ifelse(kinship > 0.1, "kin", "nonkin")) %>%
ggplot(aes(x = sri, y = resid_response)) + facet_wrap(~dyad.sex,
scales = "free_y") + geom_point(aes(color = kin), size = 2, alpha = 0.6) +
geom_smooth(method = "lm") + xlab("association rate (simple ratio index)") +
ylab("residual response call variation") + scale_colour_manual(values = c("darkred",
"grey")) + theme_bw() + theme(legend.position = "none", legend.title = element_blank())
# combine plots
(by.sex.plot <- t1 + t2 + t3 + plot_layout(ncol = 1) + plot_annotation(tag_levels = "A"))
3.3.1 Effect of kinship and association on response calling
Code
# Overdispersion test
dispersion ratio = 72.342
Pearson's Chi-Squared = 22570.653
p-value = < 0.001Code
# Overdispersion test
dispersion ratio = 0.335
Pearson's Chi-Squared = 104.319
p-value = 1# Check for zero-inflation
Observed zeros: 133
Predicted zeros: 116
Ratio: 0.87[1] 2699.808[1] 2733.723Code
[1] 2648.867[1] 2686.55 Family: nbinom2 ( log )
Formula:
n_response ~ kinship * sri + log_inquiry + offset(log(flight_time)) +
(1 | bat_flying) + (1 | bat_roosting) + (1 | group)
Zero inflation: ~1
Data: d
AIC BIC logLik deviance df.resid
2648.9 2686.6 -1314.4 2628.9 310
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
bat_flying (Intercept) 3.792e-08 0.0001947
bat_roosting (Intercept) 2.296e-01 0.4791231
group (Intercept) 1.678e-01 0.4096551
Number of obs: 320, groups: bat_flying, 50; bat_roosting, 50; group, 21
Dispersion parameter for nbinom2 family (): 0.707
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.2182 0.6294 -1.935 0.0529 .
kinship -3.3360 1.9159 -1.741 0.0817 .
sri -0.7478 0.6424 -1.164 0.2443
log_inquiry 0.1897 0.1154 1.644 0.1002
kinship:sri 5.2163 2.7487 1.898 0.0577 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Zero-inflation model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.4362 0.1242 -3.512 0.000444 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Code
[1] 2650.47[1] 2684.385 Family: nbinom2 ( log )
Formula:
n_response ~ kinship + sri + log_inquiry + offset(log(flight_time)) +
(1 | bat_flying) + (1 | bat_roosting) + (1 | group)
Zero inflation: ~1
Data: d
AIC BIC logLik deviance df.resid
2650.5 2684.4 -1316.2 2632.5 311
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
bat_flying (Intercept) 4.139e-08 0.0002035
bat_roosting (Intercept) 1.976e-01 0.4445021
group (Intercept) 1.687e-01 0.4107856
Number of obs: 320, groups: bat_flying, 50; bat_roosting, 50; group, 21
Dispersion parameter for nbinom2 family (): 0.678
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.6725 0.5816 -2.876 0.00403 **
kinship 0.2071 0.5173 0.400 0.68894
sri -0.2454 0.5926 -0.414 0.67882
log_inquiry 0.2338 0.1138 2.054 0.03995 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Zero-inflation model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.4456 0.1256 -3.548 0.000388 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Code
# get predicted values
d$n_response_predicted <- predict(fit, type = "response", newdata = d,
allow.new.levels = T)
# plot model performance
d %>%
filter(!is.na(kinship)) %>%
mutate(kinship = kinship > 0) %>%
ggplot(aes(x = n_response, y = n_response_predicted)) + geom_point(size = 2,
aes(color = kinship)) + geom_smooth(method = "lm") + xlab("observed count of response calls") +
ylab("predicted count of response calls") + scale_color_manual(values = c("darkgrey",
"darkred")) + theme(legend.position = "top")
Code
Call:
lm(formula = n_response_predicted ~ n_response, data = d)
Residuals:
Min 1Q Median 3Q Max
-105.210 -18.559 -1.269 18.411 109.081
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 58.98540 2.02304 29.16 <2e-16 ***
n_response 0.18176 0.01096 16.58 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 31.09 on 318 degrees of freedom
(126 observations deleted due to missingness)
Multiple R-squared: 0.4638, Adjusted R-squared: 0.4621
F-statistic: 275 on 1 and 318 DF, p-value: < 2.2e-16Code
# get model coefficients with 95% CIs
coefs <- tidy(fit, conf.int = TRUE, exponentiate = F, effects = "fixed",
conf.method = "profile") %>%
filter(component == "cond") %>%
dplyr::select(-effect, -component) %>%
mutate(type = "full model")
# coefs
# save model results write.csv(coefs,
# file='response_model_results.csv')
# effect of SRI with kin only
fitk <- glmmTMB(n_response ~ sri + log_inquiry + offset(log(flight_time)) +
(1 | bat_flying) + (1 | bat_roosting) + (1 | group), data = d %>%
filter(kinship > 0), ziformula = ~1, family = nbinom2)
coefs.k <- tidy(fitk, conf.int = TRUE, exponentiate = F, effects = "fixed",
conf.method = "profile") %>%
filter(component == "cond") %>%
dplyr::select(-effect, -component) %>%
mutate(type = "kin only")
# coefs.k
# effect of SRI with nonkin only
fitn <- glmmTMB(n_response ~ sri + log_inquiry + offset(log(flight_time)) +
(1 | bat_flying) + (1 | bat_roosting) + (1 | group), data = d %>%
filter(kinship == 0), ziformula = ~1, family = nbinom2)
# get 95% CI
coefs.n <- tidy(fitn, conf.int = TRUE, exponentiate = F, effects = "fixed",
conf.method = "profile") %>%
filter(component == "cond") %>%
dplyr::select(-effect, -component) %>%
mutate(type = "nonkin only")
# coefs.n
# plot model results
theme_set(theme_bw(base_size = 12))
plot1 <- coefs %>%
filter(term != "(Intercept)") %>%
mutate(term = case_when(term == "sri" ~ "association", term ==
"kinship" ~ "kinship", term == "log_inquiry" ~ "log count of inquiry calls",
term == "kinship:sri" ~ "association x kinship interaction")) %>%
mutate(term = factor(term, levels = c("log count of inquiry calls",
"association x kinship interaction", "association", "kinship"))) %>%
mutate(term = fct_rev(term)) %>%
ggplot(aes(x = estimate, y = term)) + geom_point(size = 2) + geom_errorbarh(aes(xmin = conf.low,
xmax = conf.high, height = 0.2), size = 1) + geom_vline(xintercept = 0,
linetype = "dashed") + ylab("") + xlab("coefficient") + coord_cartesian(xlim = c(-1.5,
1.5)) + theme(axis.text = element_text(size = 12), strip.text = element_text(size = 12))
plot1
Code
# save as PDF ggsave( 'response_model_results.pdf', plot =
# plot1, scale = 1, width = 5, height = 2, units = 'in', dpi =
# 300)
# check that sri and kinship do not predict response calls as
# single predictors fit zero-inflated negative binomial model
# with only sri
t <- glmmTMB(n_response ~ sri + log_inquiry + offset(log(flight_time)) +
(1 | bat_flying) + (1 | bat_roosting) + (1 | group), data = d,
ziformula = ~1, family = nbinom2)
summary(t) Family: nbinom2 ( log )
Formula:
n_response ~ sri + log_inquiry + offset(log(flight_time)) + (1 |
bat_flying) + (1 | bat_roosting) + (1 | group)
Zero inflation: ~1
Data: d
AIC BIC logLik deviance df.resid
3579.0 3611.6 -1781.5 3563.0 428
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
bat_flying (Intercept) 2.229e-08 0.0001493
bat_roosting (Intercept) 6.142e-01 0.7836871
group (Intercept) 5.587e-02 0.2363637
Number of obs: 436, groups: bat_flying, 72; bat_roosting, 73; group, 23
Dispersion parameter for nbinom2 family (): 0.646
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.1155 0.5679 -3.725 0.000195 ***
sri -0.1597 0.4814 -0.332 0.740114
log_inquiry 0.3257 0.1102 2.955 0.003122 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Zero-inflation model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.4676 0.1155 -4.048 5.17e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Code
Family: nbinom2 ( log )
Formula:
n_response ~ kinship + log_inquiry + offset(log(flight_time)) +
(1 | bat_flying) + (1 | bat_roosting) + (1 | group)
Zero inflation: ~1
Data: d
AIC BIC logLik deviance df.resid
2648.6 2678.8 -1316.3 2632.6 312
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
bat_flying (Intercept) 2.750e-08 0.0001658
bat_roosting (Intercept) 2.023e-01 0.4498245
group (Intercept) 1.786e-01 0.4226129
Number of obs: 320, groups: bat_flying, 50; bat_roosting, 50; group, 21
Dispersion parameter for nbinom2 family (): 0.681
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.7899 0.5084 -3.521 0.00043 ***
kinship 0.1260 0.4792 0.263 0.79257
log_inquiry 0.2263 0.1125 2.011 0.04429 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Zero-inflation model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.4447 0.1254 -3.545 0.000393 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Code
######## Effect of kinship and association on inquiry
######## calling---------
# poisson model is overdispersed
t <- glmer(n_inquiry ~ kinship * sri + log(n_response + 1) + offset(log(flight_time)) +
(1 | bat_flying) + (1 | bat_roosting) + (1 | group), family = poisson,
data = d)
check_overdispersion(t)# Overdispersion test
dispersion ratio = 6.731
Pearson's Chi-Squared = 2100.164
p-value = < 0.001Code
# Overdispersion test
dispersion ratio = 0.820
Pearson's Chi-Squared = 255.132
p-value = 0.991[1] 3129.521[1] 3163.436Code
# Overdispersion test
dispersion ratio = 0.813
Pearson's Chi-Squared = 252.950
p-value = 0.993[1] 3130.339[1] 3164.254 Family: nbinom2 ( log )
Formula:
n_inquiry ~ kinship * sri + n_response + offset(log(flight_time)) +
(1 | bat_flying) + (1 | bat_roosting) + (1 | group)
Data: d
AIC BIC logLik deviance df.resid
3130.3 3164.3 -1556.2 3112.3 311
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
bat_flying (Intercept) 1.833e-01 4.282e-01
bat_roosting (Intercept) 2.748e-09 5.242e-05
group (Intercept) 3.400e-02 1.844e-01
Number of obs: 320, groups: bat_flying, 50; bat_roosting, 50; group, 21
Dispersion parameter for nbinom2 family (): 4.68
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.3629986 0.1517410 -8.982 <2e-16 ***
kinship -0.5338331 0.6126168 -0.871 0.384
sri -0.1397377 0.1970890 -0.709 0.478
n_response 0.0001023 0.0002075 0.493 0.622
kinship:sri 0.4310938 0.9048669 0.476 0.634
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Code
[1] 3128.567[1] 3158.713 Family: nbinom2 ( log )
Formula:
n_inquiry ~ kinship + sri + n_response + offset(log(flight_time)) +
(1 | bat_flying) + (1 | bat_roosting) + (1 | group)
Data: d
AIC BIC logLik deviance df.resid
3128.6 3158.7 -1556.3 3112.6 312
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
bat_flying (Intercept) 1.861e-01 0.4314410
bat_roosting (Intercept) 2.571e-09 0.0000507
group (Intercept) 3.441e-02 0.1855110
Number of obs: 320, groups: bat_flying, 50; bat_roosting, 50; group, 21
Dispersion parameter for nbinom2 family (): 4.69
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.3851696 0.1446519 -9.576 <2e-16 ***
kinship -0.2524912 0.1637266 -1.542 0.123
sri -0.1002280 0.1789141 -0.560 0.575
n_response 0.0001123 0.0002065 0.544 0.587
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Code
# get predicted values
d$n_inquiry_predicted <- predict(fit2, type = "response", newdata = d,
allow.new.levels = T)
# plot model performance
d %>%
filter(!is.na(kinship)) %>%
mutate(kinship = kinship > 0) %>%
ggplot(aes(x = n_inquiry, y = n_inquiry_predicted)) + geom_point(size = 2,
aes(color = kinship)) + geom_smooth(method = "lm") + xlab("observed count of inquiry calls") +
ylab("predicted count of inquiry calls") + scale_color_manual(values = c("darkgrey",
"darkred")) + theme(legend.position = "top")
Code
Call:
lm(formula = n_inquiry_predicted ~ n_inquiry, data = d)
Residuals:
Min 1Q Median 3Q Max
-92.258 -12.967 -0.546 11.381 87.999
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 19.40700 2.06189 9.412 <2e-16 ***
n_inquiry 0.71686 0.02554 28.063 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 19.59 on 318 degrees of freedom
(126 observations deleted due to missingness)
Multiple R-squared: 0.7124, Adjusted R-squared: 0.7114
F-statistic: 787.5 on 1 and 318 DF, p-value: < 2.2e-16 Family: nbinom2 ( log )
Formula:
n_inquiry ~ kinship + sri + n_response + offset(log(flight_time)) +
(1 | bat_flying) + (1 | bat_roosting) + (1 | group)
Data: d
AIC BIC logLik deviance df.resid
3128.6 3158.7 -1556.3 3112.6 312
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
bat_flying (Intercept) 1.861e-01 0.4314410
bat_roosting (Intercept) 2.571e-09 0.0000507
group (Intercept) 3.441e-02 0.1855110
Number of obs: 320, groups: bat_flying, 50; bat_roosting, 50; group, 21
Dispersion parameter for nbinom2 family (): 4.69
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.3851696 0.1446519 -9.576 <2e-16 ***
kinship -0.2524912 0.1637266 -1.542 0.123
sri -0.1002280 0.1789141 -0.560 0.575
n_response 0.0001123 0.0002065 0.544 0.587
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Code
coefs2 <- tidy(fit2, conf.int = TRUE, exponentiate = F, effects = "fixed",
conf.method = "profile") %>%
filter(component == "cond") %>%
dplyr::select(-effect, -component) %>%
mutate(type = "full model")
# coefs2
# save model results write.csv(coefs2,
# file='inquiry_model_results.csv')
# plot model results
theme_set(theme_bw(base_size = 12))
plot2 <- coefs2 %>%
filter(term != "(Intercept)") %>%
mutate(term = case_when(term == "sri" ~ "association", term ==
"kinship" ~ "kinship", term == "n_response" ~ "response calls")) %>%
mutate(term = factor(term, levels = c("response calls", "association",
"kinship"))) %>%
mutate(term = fct_rev(term)) %>%
ggplot(aes(x = estimate, y = term)) + geom_point(size = 2) + geom_errorbarh(aes(xmin = conf.low,
xmax = conf.high, height = 0.2), size = 1) + geom_vline(xintercept = 0,
linetype = "dashed") + ylab("") + xlab("estimate (log odds)") +
theme(axis.text = element_text(size = 12), strip.text = element_text(size = 12))
plot2
Code
# save as PDF ggsave( 'inquiry_model_results.pdf', plot = plot1,
# scale = 1, width = 7, height = 3.5, units = 'in', dpi = 300)
# get relative amount of variance explained by random intercept
tidy(fit2, conf.int = F, exponentiate = F) %>%
filter(effect == "ran_pars") %>%
mutate(variance = estimate^2) %>%
dplyr::select(effect, group, variance) %>%
group_by(effect) %>%
mutate(sum = sum(variance)) %>%
mutate(prop = variance/sum)# A tibble: 3 × 5
# Groups: effect [1]
effect group variance sum prop
<chr> <chr> <dbl> <dbl> <dbl>
1 ran_pars bat_flying 0.186 0.221 0.844
2 ran_pars bat_roosting 0.00000000257 0.221 0.0000000117
3 ran_pars group 0.0344 0.221 0.156 Code
Family: nbinom2 ( log )
Formula:
n_inquiry ~ sri + n_response + offset(log(flight_time)) + (1 |
bat_flying) + (1 | bat_roosting) + (1 | group)
Data: d
AIC BIC logLik deviance df.resid
4257.0 4285.5 -2121.5 4243.0 429
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
bat_flying (Intercept) 0.23098 0.48060
bat_roosting (Intercept) 0.00448 0.06694
group (Intercept) 0.01246 0.11164
Number of obs: 436, groups: bat_flying, 72; bat_roosting, 73; group, 23
Dispersion parameter for nbinom2 family (): 5.1
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.3782473 0.1046922 -13.165 <2e-16 ***
sri -0.1730846 0.1317072 -1.314 0.1888
n_response 0.0003027 0.0001691 1.790 0.0734 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Code
Family: nbinom2 ( log )
Formula:
n_inquiry ~ kinship + n_response + offset(log(flight_time)) +
(1 | bat_flying) + (1 | bat_roosting) + (1 | group)
Data: d
AIC BIC logLik deviance df.resid
3126.9 3153.3 -1556.4 3112.9 313
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
bat_flying (Intercept) 1.825e-01 4.272e-01
bat_roosting (Intercept) 3.244e-09 5.695e-05
group (Intercept) 3.491e-02 1.868e-01
Number of obs: 320, groups: bat_flying, 50; bat_roosting, 50; group, 21
Dispersion parameter for nbinom2 family (): 4.67
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.4480670 0.0905664 -15.989 <2e-16 ***
kinship -0.2738852 0.1594311 -1.718 0.0858 .
n_response 0.0001283 0.0002052 0.625 0.5319
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Takeaways
- Individual ID but not kinship or association explain calling rates
Session information
R version 4.2.2 Patched (2022-11-10 r83330)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.5 LTS
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
locale:
[1] LC_CTYPE=es_ES.UTF-8 LC_NUMERIC=C
[3] LC_TIME=es_CR.UTF-8 LC_COLLATE=es_ES.UTF-8
[5] LC_MONETARY=es_CR.UTF-8 LC_MESSAGES=es_ES.UTF-8
[7] LC_PAPER=es_CR.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] MASS_7.3-58.2 performance_0.10.3 patchwork_1.1.2
[4] boot_1.3-28 broom.mixed_0.2.9.4 glmmTMB_1.1.7
[7] lme4_1.1-33 Matrix_1.5-1 mapproj_1.2.11
[10] maps_3.4.1 geosphere_1.5-18 lubridate_1.9.2
[13] forcats_1.0.0 stringr_1.5.0 dplyr_1.1.0
[16] purrr_1.0.1 readr_2.1.4 tidyr_1.3.0
[19] tibble_3.2.1 tidyverse_2.0.0 RColorBrewer_1.1-3
[22] statnet_2019.6 tsna_0.3.5 ergm.count_4.1.1
[25] tergm_4.2.0 networkDynamic_0.11.3 ergm_4.5.0
[28] igraphdata_1.0.1 ecodist_2.0.9 ggmap_3.0.1
[31] ggplot2_3.4.3 assortnet_0.20 asnipe_1.1.17
[34] bipartite_2.18 sna_2.7-1 network_1.18.1
[37] statnet.common_4.9.0 vegan_2.6-4 lattice_0.20-45
[40] permute_0.9-7 igraph_1.5.1
loaded via a namespace (and not attached):
[1] minqa_1.2.5 colorspace_2.1-0 ellipsis_0.3.2
[4] rprojroot_2.0.3 estimability_1.4.1 xaringanExtra_0.7.0
[7] rstudioapi_0.14 farver_2.1.1 listenv_0.8.0
[10] furrr_0.3.1 remotes_2.4.2 mvtnorm_1.1-3
[13] fansi_1.0.4 codetools_0.2-19 splines_4.2.2
[16] cachem_1.0.8 robustbase_0.99-0 knitr_1.43
[19] spam_2.9-1 jsonlite_1.8.7 nloptr_2.0.3
[22] packrat_0.8.1 broom_1.0.4 cluster_2.1.4
[25] png_0.1-8 compiler_4.2.2 httr_1.4.6
[28] backports_1.4.1 emmeans_1.8.6 fastmap_1.1.1
[31] cli_3.6.1 formatR_1.12 htmltools_0.5.6
[34] tools_4.2.2 dotCall64_1.0-2 coda_0.19-4
[37] gtable_0.3.4 glue_1.6.2 Rcpp_1.0.11
[40] rle_0.9.2 vctrs_0.6.3 nlme_3.1-162
[43] insight_0.19.2 xfun_0.40 globals_0.16.1
[46] rbibutils_2.2.15 trust_0.1-8 timechange_0.2.0
[49] lifecycle_1.0.3 future_1.28.0 DEoptimR_1.1-2
[52] scales_1.2.1 hms_1.1.2 parallel_4.2.2
[55] ergm.multi_0.2.0 TMB_1.9.6 fields_15.2
[58] yaml_2.3.7 memoise_2.0.1 stringi_1.7.12
[61] RgoogleMaps_1.4.5.3 Rdpack_2.5 rlang_1.1.1
[64] pkgconfig_2.0.3 bitops_1.0-7 evaluate_0.21
[67] lpSolveAPI_5.5.2.0-17.10 labeling_0.4.3 htmlwidgets_1.5.4
[70] tidyselect_1.2.0 parallelly_1.32.1 plyr_1.8.7
[73] magrittr_2.0.3 R6_2.5.1 generics_0.1.3
[76] pillar_1.9.0 withr_2.5.0 mgcv_1.8-41
[79] sp_2.0-0 crayon_1.5.2 utf8_1.2.3
[82] networkLite_1.0.5 tzdb_0.3.0 rmarkdown_2.21
[85] jpeg_0.1-9 grid_4.2.2 sketchy_1.0.2
[88] digest_0.6.33 xtable_1.8-4 numDeriv_2016.8-1.1
[91] munsell_0.5.0 viridisLite_0.4.2